function F = cal_model_moment(param, const, eqm)

    % Read-in
    sigma           = param.sigma;
    d               = param.d;
    Q               = param.Q;    
    gamma           = const.gamma;
    phi_v           = const.phi_v;
    phi_y           = const.phi_y;
    lnzQ            = const.lnzQ;        
    density_omega   = const.density_omega;    
    lnwage_grid     = const.lnwage_grid;  
    q_indicator     = eqm.q_indicator;
    nq              = eqm.nq;
    nd              = eqm.nd;
    quintileflag    = eqm.quintileflag;   
    E_zx00           = eqm.integ.E_zx00;
    E_zx01           = eqm.integ.E_zx01;
    E_zx10           = eqm.integ.E_zx10;
    E_zx11           = eqm.integ.E_zx11;           
    ExportProbQ     = eqm.ExportProbQ;   
    ImportProbQ1    = eqm.ImportProbQ1;
    ImportProbQ0    = eqm.ImportProbQ0;
    Pi00            = eqm.Pi00;
    Pi01            = eqm.Pi01;
    Pi10            = eqm.Pi10;
    Pi11            = eqm.Pi11;
    V               = eqm.V;
    M               = eqm.M;
    theta_v         = eqm.theta_v;
    theta_m         = eqm.theta_m; 
    X_m             = eqm.X_m;
    Psig            = eqm.Psig;
    c_H             = eqm.c_H;
    rv_rev          = eqm.rv_rev;   
    rm_spend        = eqm.rm_spend;
    Vbar            = eqm.Vbar;
    

%% Useful Notation

    % n(q,E)
    nq1 = density_omega'*(ExportProbQ.*q_indicator);
    nq0 = density_omega'*((1-ExportProbQ).*q_indicator);

    % n(q,E,I)
    nq11 = density_omega'*(ExportProbQ.*ImportProbQ1.*q_indicator);
    nq10 = density_omega'*(ExportProbQ.*(1-ImportProbQ1).*q_indicator);
    nq01 = density_omega'*((1-ExportProbQ).*ImportProbQ0.*q_indicator);
    nq00 = density_omega'*((1-ExportProbQ).*(1-ImportProbQ0).*q_indicator);

    % E_lnz(q,E,I)
    E_lnz11 = density_omega'*(lnzQ.*ExportProbQ.*ImportProbQ1.*q_indicator);
    E_lnz10 = density_omega'*(lnzQ.*ExportProbQ.*(1-ImportProbQ1).*q_indicator);
    E_lnz01 = density_omega'*(lnzQ.*(1-ExportProbQ).*ImportProbQ0.*q_indicator);
    E_lnz00 = density_omega'*(lnzQ.*(1-ExportProbQ).*(1-ImportProbQ0).*q_indicator);

    % E_[lnz^2](q,E,I)
    E2_lnz11 = density_omega'*(lnzQ.^2.*ExportProbQ.*ImportProbQ1.*q_indicator);
    E2_lnz10 = density_omega'*(lnzQ.^2.*ExportProbQ.*(1-ImportProbQ1).*q_indicator);
    E2_lnz01 = density_omega'*(lnzQ.^2.*(1-ExportProbQ).*ImportProbQ0.*q_indicator);
    E2_lnz00 = density_omega'*(lnzQ.^2.*(1-ExportProbQ).*(1-ImportProbQ0).*q_indicator);


%% Mean In-degree and Out-degree

    % Mean degree by q
    mean_indegree_q = M.*theta_m./nq;    
    mean_outdegree_q = Vbar.*sum(phi_v.*repmat(theta_v',1,Q),1)./nq;   
    mean_indegree_q(nq==0)=0;    
    mean_outdegree_q(nq==0)=0;     

    % Mean degree by quintile
    mean_indegree_d =  (sum(repmat(mean_indegree_q.*nq,d,1).*quintileflag,2)./sum(repmat(nq,d,1).*quintileflag,2))';    
    mean_outdegree_d =(sum(repmat(mean_outdegree_q.*nq,d,1).*quintileflag,2)./sum(repmat(nq,d,1).*quintileflag,2))';  
         
    % Mean degrees for all
    mean_indegree = sum(mean_indegree_q.*nq);   
    mean_outdegree = sum(mean_outdegree_q.*nq);
       
    % Store
    model.mean_indegree_q           = mean_indegree_q;
    model.mean_indegree_d           = mean_indegree_d;
    model.mean_indegree             = mean_indegree;  
    model.mean_outdegree_q          = mean_outdegree_q;
    model.mean_outdegree_d          = mean_outdegree_d;
    model.mean_outdegree            = mean_outdegree;


%% Mean and Variance of Log Sales

    % Mean of log sales
    mean_lnsales_q = ( (sigma-1)*gamma*E_lnz11+log(Pi11).*nq11 ...
                      +(sigma-1)*gamma*E_lnz10+log(Pi10).*nq10 ...
                      +(sigma-1)*gamma*E_lnz01+log(Pi01).*nq01 ...
                      +(sigma-1)*gamma*E_lnz00+log(Pi00).*nq00 )./nq;                  
    mean_lnsales_q(nq==0)=0;   
    mean_lnsales_d = (sum(repmat(mean_lnsales_q.*nq,d,1).*quintileflag,2)./sum(repmat(nq,d,1).*quintileflag,2))';
    mean_lnsales = sum(mean_lnsales_q.*nq);

    % Variance of log sales
    var_lnsales_q =( ((sigma-1)*gamma)^2*E2_lnz11+2*(sigma-1)*gamma*log(Pi11).*E_lnz11+(log(Pi11)).^2.*nq11 ...
                    +((sigma-1)*gamma)^2*E2_lnz10+2*(sigma-1)*gamma*log(Pi10).*E_lnz10+(log(Pi10)).^2.*nq10 ...
                    +((sigma-1)*gamma)^2*E2_lnz01+2*(sigma-1)*gamma*log(Pi01).*E_lnz01+(log(Pi01)).^2.*nq01 ...
                    +((sigma-1)*gamma)^2*E2_lnz00+2*(sigma-1)*gamma*log(Pi00).*E_lnz00+(log(Pi00)).^2.*nq00 ...
                   )./nq - (mean_lnsales_q).^2;
    var_lnsales_q(nq==0)=0;   
    var_lnsales_d = (sum(repmat((var_lnsales_q+mean_lnsales_q.^2).*nq,d,1).*quintileflag,2)./sum(repmat(nq,d,1).*quintileflag,2))'-mean_lnsales_d.^2;   
    var_lnsales = sum((var_lnsales_q+mean_lnsales_q.^2).*nq)-mean_lnsales.^2;  

    % Store
    model.mean_lnsales_q = mean_lnsales_q;
    model.mean_lnsales_d = mean_lnsales_d;
    model.mean_lnsales = mean_lnsales;
    model.var_lnsales_q = var_lnsales_q;
    model.var_lnsales_d = var_lnsales_d;
    model.var_lnsales = var_lnsales; 


%% Share of Total Network Sales
   
    % Mean network sales
    mean_network_sales_q = (Pi00.*E_zx00+rv_rev.*Pi10.*E_zx10 ...
                           +Pi01.*E_zx01+rv_rev.*Pi11.*E_zx11)./nq;
    mean_network_sales_q(nq==0)=0;
    mean_network_sales_d = (sum(repmat(mean_network_sales_q.*nq,d,1).*quintileflag,2)./sum(repmat(nq,d,1).*quintileflag,2))'; 

    % Network sales share
    network_sales_share_d = mean_network_sales_d.*nd/sum(mean_network_sales_d.*nd);        
        
    % Store 
    model.mean_network_sales_q  = mean_network_sales_q;
    model.mean_network_sales_d  = mean_network_sales_d;
    model.network_sales_share_d = network_sales_share_d;    
        

%% Exporter and Importer Share

    % Fraction of exporters 
    exporter_share_q = nq1./nq;
    exporter_share_q(nq==0)=0;
    exporter_share_d = (sum(repmat(exporter_share_q.*nq,d,1).*quintileflag,2)./sum(repmat(nq,d,1).*quintileflag,2))';
    exporter_share = sum(exporter_share_q.*nq);  

    % Fraction of importers
    importer_share_q = (nq11+nq01)./nq;
    importer_share_q(nq==0)=0;
    importer_share_d = (sum(repmat(importer_share_q.*nq,d,1).*quintileflag,2)./sum(repmat(nq,d,1).*quintileflag,2))';
    importer_share = sum(importer_share_q.*nq);         

    % Store
    model.exporter_share_q      = exporter_share_q;
    model.exporter_share_d      = exporter_share_d;
    model.exporter_share        = exporter_share;
    model.importer_share_q      = importer_share_q;
    model.importer_share_d      = importer_share_d;
    model.importer_share        = importer_share;
    
    
%% Export and Import Intensity

    % Export intensity of exporters
    export_intensity_q = 1-rv_rev;
    export_intensity_d = (sum(repmat((1-rv_rev).*exporter_share_q.*nq,d,1).*quintileflag,2)./sum(repmat(exporter_share_q.*nq,d,1).*quintileflag,2))';
    export_intensity = sum((1-rv_rev).*nq); 

    % Import intensity of importers
    import_intensity_q = 1-rm_spend;
    import_intensity_d = (sum(repmat((1-rm_spend).*importer_share_q.*nq,d,1).*quintileflag,2)./sum(repmat(importer_share_q.*nq,d,1).*quintileflag,2))';
    import_intensity = sum((1-rm_spend).*nq);

    % Store
    model.export_intensity_q    = export_intensity_q;
    model.export_intensity_d    = export_intensity_d;
    model.export_intensity      = export_intensity;
    model.import_intensity_q    = import_intensity_q;
    model.import_intensity_d    = import_intensity_d;
    model.import_intensity      = import_intensity;

    
%% Average Log Wage of Suppliers

    % Unweighted average log wage of suppliers
    avg_unwgt_lnwageS_q = 1./V.*sum(phi_v.*repmat((log(param.w)+lnwage_grid).*Vbar, Q,1), 2)';    
    avg_unwgt_lnwageS_d = (sum(repmat(avg_unwgt_lnwageS_q.*nq, d, 1).*quintileflag, 2)./sum(repmat(nq,d,1).*quintileflag,2))';   

    % Weighted average log wage of suppliers
    avg_wgt_lnwageS_q = theta_m.*c_H.^(sigma-1)./V.*sum(phi_y.*phi_v.*repmat((log(param.w)+lnwage_grid).*Psig, Q,1),2)';  
    avg_wgt_lnwageS_d = (sum(repmat(avg_wgt_lnwageS_q.*nq, d, 1).*quintileflag, 2)./sum(repmat(nq,d,1).*quintileflag,2))';  
    
    % Store    
    model.avg_wgt_lnwageS_q     = avg_wgt_lnwageS_q;   
    model.avg_wgt_lnwageS_d     = avg_wgt_lnwageS_d;
    model.avg_wgt_lnwageS_d4    = avg_wgt_lnwageS_d(2:5)-avg_wgt_lnwageS_d(1);    
    model.avg_unwgt_lnwageS_q   = avg_unwgt_lnwageS_q;
    model.avg_unwgt_lnwageS_d   = avg_unwgt_lnwageS_d;
    model.avg_unwgt_lnwageS_d4  = avg_unwgt_lnwageS_d(2:5)-avg_unwgt_lnwageS_d(1);
    

%% Untargeted: Aggregate Network Moments

    % Supplier and customer shares
    common = mysum(repmat(theta_v',1,Q).*phi_v.*repmat(Vbar,Q,1), quintileflag, param);
    share_indegree_d = common./repmat(sum(common,2),1,d);
    share_outdegree_d = common./repmat(sum(common,1),d,1);

    % Expenditure and sales share
    common = mysum(repmat((theta_m./V.*c_H.^(sigma-1).*X_m)',1,Q).*phi_y.*phi_v.*repmat(Psig,Q,1), quintileflag, param);
    share_spend_d = common./repmat(sum(common,2),1,d);
    share_sales_d = common./repmat(sum(common,1),d,1);
    
    % Store
    model.share_indegree_d  = share_indegree_d;
    model.share_outdegree_d = share_outdegree_d;
    model.share_spend_d     = share_spend_d;
    model.share_sales_d     = share_sales_d;
    
   
%% Untargeted: Wage Regression Coefficients

    % Decompose weighted average log wage of suppliers
    avg_wgt_lnwageS_ext_q = avg_unwgt_lnwageS_q;
    avg_wgt_lnwageS_int_q = avg_wgt_lnwageS_q - avg_unwgt_lnwageS_q;    
    
    % Wage regression
    reg_W = diag(nq);
    reg_x = [ones(Q,1), param.w.*lnwage_grid'];
    reg_y_ext = avg_wgt_lnwageS_ext_q';
    reg_y_int = avg_wgt_lnwageS_int_q';
    reg_y_sum = avg_wgt_lnwageS_q';
    beta_ext = (reg_x'*reg_W*reg_x)\(reg_x'*reg_W*reg_y_ext);
    beta_int = (reg_x'*reg_W*reg_x)\(reg_x'*reg_W*reg_y_int);
    beta_sum = (reg_x'*reg_W*reg_x)\(reg_x'*reg_W*reg_y_sum);
    slope_ext = beta_ext(2);
    slope_int = beta_int(2);
    slope_sum = beta_sum(2);   

    % Store    
    model.avg_wgt_lnwageS_ext_q = avg_wgt_lnwageS_ext_q;
    model.avg_wgt_lnwageS_int_q = avg_wgt_lnwageS_int_q;
    model.slope_ext             = slope_ext;
    model.slope_int             = slope_int;
    model.slope_sum             = slope_sum;    
  

%% RETURN
    
    F = model;

    
end